/**
* @Arquivo: lagrangian.c
* 
*	Geração e Solução do problema de coloração usando a formulção de representantes 
*
* @author: Marcio Costa Santos
* @data: 20/05/2010
*
*/

#include"lagrangian.h"
#include<stdio.h>
#include<stdlib.h>
#include <string.h>


double weigths[MAX_VALUE][MAX_VALUE];/*o vetor de pesos*/

int status;/*o status do problema*/

extern int n;/*o numero de vertices*/

extern int m;/*o numero de arestas*/

extern long va[];/*a ordem dos verties*/

extern FILE* out;/*o arquivo de saida*/

long edgeVar[MAX_VALUE][MAX_VALUE];/*correspondencia verticeXvertice -> aresta*/

long vertexVar[MAX_VALUE][MAX_VALUE];/* correspondencia aresta -> vertexXvertex*/

long var[MAX_VALUE] ;/*numero de variaveis*/

int sink[MAX_VALUE];

double subpobj;

CPXENVptr env[MAX_VALUE];/*o ambiente do problema*/

CPXLPptr lp[MAX_VALUE];/*os problemas*/

double x[MAX_VALUE][MAX_VALUE];/*o vetor de pesos*/

double bestx[MAX_VALUE][MAX_VALUE];

double bestLB;/*o melhor limite inferior encontrado*/

double bestUB;/*o melhor limite superior encontrado*/

double y[MAX_VALUE];/*vetor auxiliar para recuperar solução*/

double aux[MAX_VALUE];/*vetor auxiliar para recuperar solução*/

double sol[MAX_VALUE][MAX_VALUE];/*uma solução viavel do problema*/

double lb;/*o limite inferior*/

double ub;/*o limite superior*/

double pi;/*passo do lagrangeano*/

int noI;/*melhoria?*/

int handle;/*auxiliar*/

double coef[MAX_VALUE];/*os coeficientes da função objetivo*/

int indices[MAX_VALUE];/*auxiliar para os coeficientes da função objetivo*/

char buffer[100];/*string para controle de erro*/

/**
*	Verifica se as rotinas do Cplex contem algum erro pelo valor retorndo
*/
int getErro(char* sucessMsg, char* failMsg, int v){
  if(status == 0){}
  else{
      if ( CPXgeterrorstring (env[v], status, buffer) != NULL ) printf ("%s \n", buffer);
      else printf ("CPLEX Error %3d:  codigo de erro desconhecido.\n",status);
      printf ("%s \n", failMsg);
  }
  return status;
}

/** 
	Inicia os vetores usados no programa
*/
int prepareVariables(){
  int i, j;
  pi = 4.0;
  long ckeck = 0;
  noI = 0;
  lb = ub =  0;
  bestLB = (-1)*n;
  bestUB = n;
  for(i = 0; i < n; i++){
      mu[i] = Lx[i] = 1.0;
      y[i] = 0.0;
      indices[i] = i;
      for(j = 0; j < n; j++){
          edgeVar[i][j] = 0;
          vertexVar[i][j] = 0;
          x[i][j] = 0.0;
      }
  }

  for(i = 0; i < n; i++){
      var[va[i]] = 0;
      ckeck = 1;
      for(j = i + 1; j < n; j++){
          if(hasEdge(va[i],va[j])){
              edgeVar[va[i]][va[j]] = var[va[i]];
              vertexVar[va[i]][var[va[i]]]= va[j];
              var[va[i]]++;
              ckeck = 0;
          }
      }
      sink[va[i]] = ckeck;
  }
  return 0;
}

/**
	Adicona as restrições para cada vértice
*/
int addConstrainsForVertex(long v){
  int i, j;
  long constrains = 0;/*numero de restrições*/
  long nonNulosCoef =0;/*numero de coeficientes não nulos*/
  double values[MAX_VALUE*MAX_VALUE];/*os valores nao nulos*/
  int valuesColumns[MAX_VALUE*MAX_VALUE];/*variavel para montar o problema no CPLEX*/
  int valuesRows[MAX_VALUE];/*variavel para montar o problema no CPLEX*/
  double b[MAX_VALUE];/*o lado direito das restricoes*/
  char sense[MAX_VALUE];/*o sinal da desigualdade*/
  for(i = v + 1; i < n; i++){
      if(hasEdge(va[v],va[i])){
          for(j = i + 1; j < n; j++){
              if((i < j)&&(!hasEdge(va[i], va[j]))&&(hasEdge(va[v], va[j]))){
                  /*x(v,w) + x(v,u) <= 1*/
                  valuesRows[constrains] = nonNulosCoef;
                  constrains++;
                  values[nonNulosCoef] = 1.0;
                  valuesColumns[nonNulosCoef] = edgeVar[va[v]][va[i]];
                  values[nonNulosCoef+1] = 1.0;
                  valuesColumns[nonNulosCoef+1] = edgeVar[va[v]][va[j]];
                  nonNulosCoef += 2;
              }
          }
      }
  }
  valuesRows[constrains] = nonNulosCoef;
  for(i = 0; i< constrains; i++){
      b[i] = 1.0;
      sense[i] = 'L';
  }
  status = CPXaddrows(env[va[v]], lp[va[v]], 0, constrains, nonNulosCoef, b, sense, valuesRows, valuesColumns, values, NULL, NULL);
  CPXwriteprob (env[va[v]], lp[va[v]], "myprob.lp", NULL);
  return getErro("Restriçoes adicionadas com sucesso","Problemas ao adicionar restrições",v);
}


/**
* Constroi o problema de programação linear
*/
int makePL(long v){
  int i;
  env[va[v]] = CPXopenCPLEX(&status);/*Inicia o ambiente*/
  if(getErro("ILOG/CPLEX iniciado com sucesso","Problemas para abrir o ILOG/CPLEX",v)==0){
      lp[va[v]] = CPXcreateprob(env[va[v]], &status,"Formulacao de Representantes");
      if(getErro("Problema criado com sucesso", "Problemas para criar o problema",v)==0){
          double lb[MAX_VALUE];/*o limite inferior das variaveis*/
          double ub[MAX_VALUE];/*o limite superior das variaveis*/
          char type[MAX_VALUE];/*o tipo das variaveis*/
          double obj[MAX_VALUE];
          for(i = 0; i < var[va[v]]; i++){
              lb[i] = 0.0;
              ub[i] = 1.0;
              obj[i] = 0.0;
              type[i] = 'B';
          }
          status = CPXnewcols (env[va[v]], lp[va[v]], var[va[v]], obj, lb, ub, type, NULL);
          if(addConstrainsForVertex(v)!= 0) return 1;
          if(getErro("Problema montado com sucesso", "Problemas para montado o problema",v)==0) return 1;
      }
   }
   return 0;
}

int updateX(long v){
  int i;
  for(i = 0; i < var[va[v]]; i++){
      x[va[v]][vertexVar[va[v]][i]] = y[i];
  }
  x[va[v]][va[v]] = 1.0;
  return 0;
}

int zeraX(long v){
  int i;
  for(i = 0; i < n; i++){
      x[va[v]][va[i]] = 0.0;
  }
  x[va[v]][va[v]] = 0.0;
  return 0;
}

int color_G(int *missed,const int miss){
  int i,j;
  int Rmissed[miss+1];
  int Rmiss = 0;
  if(miss == 0) return 0;
  else{
      sol[missed[0]][missed[0]] = 1;
      for(i = 1; i< miss; i++ ){
          if(missed[i] > -1){
              if(hasEdge(missed[0], missed[i])){
                  sol[missed[0]][missed[i]] = 1;
                  for(j = i; j < miss; j++){
                      if(missed[j] > -1){
                          if(!hasEdge(missed[i],missed[j])){
                              Rmissed[Rmiss] = missed[j];
                              Rmiss++;
                              missed[j] = -1;
                          }
                      }
                  }
              }else{
                  Rmissed[Rmiss] = missed[i];
                  Rmiss++;
                  missed[i] = -1;
              }
          }
      }
  }
  color_G(Rmissed, Rmiss);
  return 0;
}


double calculateUB(){
  int i,j,k;
  double value = 0.0;
  for(i = 0; i < n; i++){
      CPXgetx(env[va[i]], lp[va[i]], y, 0, CPXgetnumcols(env[va[i]], lp[va[i]])-1);
      for(j = 0; j < n; j++){
          sol[va[i]][va[j]] = 0.0;
      }
      sol[va[i]][va[i]] = x[va[i]][va[i]];
      if(sol[va[i]][va[i]] > 0.0){
          for(j = 0; j < var[va[i]]; j++){
              sol[va[i]][vertexVar[va[i]][j]] = y[j];
          }
      }
  }
  double verif;
  int miss = 0;
  int missed[MAX_VALUE];
  for(i = 0; i < n; i++){
      verif = 0.0;
      for(j = 0; j <= i; j++){
          verif += sol[va[j]][va[i]];
      }
      if(verif < 0.9){
          missed[miss] = va[i];
          miss++;
      }
  }
  color_G(missed, miss);
  for(i = 0; i < n; i++){
      value += sol[va[i]][va[i]];
  }
  if(value < bestUB){
      bestUB = value;
  }
  return value;
}

int noImproviment(){
  noI++;
  if(noI > 20){
      pi = pi/2.0;
      noI = 0;
  }
  return 0;
}


double calculateLB(){
  int i,j;
  double value = 0.0;
  for(i = 0; i < n; i++){
      value += Lx[va[i]];
  }
  for(i = 0; i < n; i++){
      for(j = i ; j < n; j++){
          if(x[va[i]][va[j]] > 0.0){
              printf("x[%d][%d] = %f\n", va[i]+1, va[j]+1, x[va[i]][va[j]]);
          }
      }
  }
  if(value > bestLB){
      printf("Atualizei o lb para %f \n", value);
      bestLB = value;
      for(i = 0; i < n; i++){
          for(j = 0 ; j < n; j++){
              bestx[va[i]][va[j]] = x[va[i]][va[j]];
          }
      }
  }
  else noImproviment();
  return value;
}

	
int updateMuLambda(){
  double g[MAX_VALUE];
  double theta = 0.0;
  double theta_div = 0.0;
  int i,j;
  for(i = 0; i < n; i++){
      g[va[i]] = 1.0;
      for(j = 0; j <= i; j++){
          g[va[i]] -= x[va[j]][va[i]];
      }
      theta_div += g[va[i]]*g[va[i]];
  }
  if(theta_div >=  0.001)theta = pi*(bestUB - lb)/theta_div;
  else theta = pi/n;
  for(i = 0; i < n; i++){
      mu[va[i]] += theta*g[va[i]];
      if(mu[va[i]] < 0.0) mu[va[i]] = 0.0;
  }
  return 0;
}

/**
*	Adiciona uma função callbaxk para poder recuperar os valores reais
*/
int mycallback (CPXCENVptr envRelax, void  *cbdata, int wherefrom, void *cbhandle, int type, int sos, int nodecnt, int bdcnt, const double  *nodeest, const int *nodebeg, const int *indices, const char *lu, const int *bd, int *useraction_p){
  int j;
  status = CPXgetcallbacknodeobjval (envRelax, cbdata, wherefrom, &subpobj);
  status = CPXgetcallbacknodex (envRelax, cbdata,wherefrom, y, 0, var[handle] -1);
  getErro("Valores do problema adquiridos com sucesso na callback","Problemas para recupara a solucao do problema na callback", handle);
  return 0;
}

int mycallback2(CPXCENVptr envRelax, void *cbdata, int wherefrom, void *cbhandle, int *useraction_p){
  int j;
  status = CPXgetcallbacknodeobjval (envRelax, cbdata, wherefrom, &subpobj);
  status = CPXgetcallbacknodex (envRelax, cbdata, wherefrom, y, 0, var[handle] -1);
  getErro("Valores do problema adquiridos com sucesso na callback","Problemas para recupara a solucao do problema na callback", handle);
  return 0;
}


int setParam(){
  int i;
  for(i = 0; i < n; i++){
      CPXsetintparam(env[i], CPX_PARAM_MIPCBREDLP, CPX_OFF);
      CPXsetintparam(env[i], CPX_PARAM_DISJCUTS , 2);
      CPXsetintparam(env[i], CPX_PARAM_FLOWCOVERS  , 2);
      CPXsetintparam(env[i], CPX_PARAM_FLOWPATHS  , 2);
      CPXsetintparam(env[i], CPX_PARAM_FRACCUTS  , 2);
      CPXsetintparam(env[i], CPX_PARAM_GUBCOVERS  , 2);
      CPXsetintparam(env[i], CPX_PARAM_IMPLBD   , 2);
      CPXsetintparam(env[i], CPX_PARAM_MIRCUTS  , 2);
      CPXsetintparam(env[i], CPX_PARAM_SUBALG , 2);
      CPXsetintparam(env[i], CPX_PARAM_CLIQUES , 2);
      CPXsetintparam(env[i], CPX_PARAM_COVERS  , 2);
      CPXsetintparam(env[i], CPX_PARAM_NODELIM, 0);/*opção para limitar a arvore de branch*/
  }
  return status;

}

/**
*	Seleciona algumas opções do cplex e resolve o problema
*/
double solvePL(int v){
  handle = va[v];
  CPXsetcutcallbackfunc(env[va[v]], mycallback2, va[v]);
  CPXsetbranchcallbackfunc(env[va[v]], mycallback, va[v]);
  CPXchgobjsen(env[va[v]], lp[va[v]], CPX_MAX);
  if(CPXmipopt(env[va[v]],lp[va[v]]) == CPXMIP_OPTIMAL){
      CPXgetx(env[va[v]],lp[va[v]], y,0, var[va[v]] -1 );
      CPXgetbestobjval(env[va[v]], lp[va[v]], &subpobj);
  }
  if(sink[va[v]] == 1 ){
      zeraX(v);
      if(mu[va[v]] < 1.0){
          x[va[v]][va[v]] = 0;
          return mu[va[v]];
      }
      else{
          x[va[v]][va[v]] = 1;
          return 1.0 ;
      }
  }
  if(1.0 - subpobj < mu[va[v]] ) {
      updateX(v);
      x[va[v]][va[v]] = 1.0;
      return 1.0 - subpobj;
  }
  else{
      zeraX(v);
      x[va[v]][va[v]] = 0.0;
      return mu[va[v]];
  }
}

int addObjFunc(long v){
  int i;
  for(i = 0; i < var[va[v]]; i++){
      coef[i] = mu[vertexVar[va[v]][i]];
  }
  status = CPXchgobj (env[va[v]], lp[va[v]], var[va[v]], indices, coef);
  CPXwriteprob (env[va[v]], lp[va[v]], "myprob.lp", NULL);
  return 0;
}


double lagrangian(double* sol){
  int i,j;
  prepareVariables();
  biuldSelectionProblem();
  for(i = 0; i < n; i++){
      makePL(i);
  }
  setParam();
  while(pi > 0.005){
      for(i = 0; i < n; i++){
          addObjFunc(i);
          Lx[va[i]] = solvePL(i);
      }
      lb = calculateLB();
      ub = calculateUB();
      updateMuLambda();
  }
  *sol = bestUB;
  return bestLB;
}

/**
*	Recupera os valores das variaveis e libera a memoria usada
*/
int getValues(){
  int valid = 0;
  int i,j;
  for(i = 0; i< n; i++){
      for(j = 0; j< n; j++){
          weigths[va[i]][va[j]] = bestx[va[i]][va[j]];
      }
      CPXfreeprob(env[i],&(lp[i]));
      CPXcloseCPLEX(&(env[i]));
  }
  return 0;
}


